#!/usr/bin/perl -w
##This script run fastq-mcf to do adapter sequence trimming (also trimmed bad base quality <= 2 as BSMAP did) on fastq files and then do statistics on adapter contaminated number.
## Allow 10% mismatches of adapter sequences, requires at least 5 bases of adapter sequences, when adapter contamination is more than 0.25%, it will begin to do trimming. 
## skewed option is disabled, quality trimming is <=2
## Total reads number is actually read1+read2 for PE.
## Adapter dimmer here is actually the reads number that is too short after clipping, but most of time it is because of adapter dimmer, not due to quality score trimming.


## author: Yaping Liu  lyping1986@gmail.com 
## time: 2013-1-9

#Usage:  perl wrap_fastq-mcf.pl [option] 1st_end.fastq [2nd_end.fastq]

my $FASTQMCF = "/export/uec-gs1/laird/users/yaping/software/ea-utils.1.1.2-484/fastq-mcf";
my $ADAPTER = "/export/uec-gs1/laird/users/yaping/data/genome_data/illumina_adapter_new.fa";


my $read1 = $ARGV[0] || die "no input reads specified";
my $read2 = $ARGV[1];

my $out1 = $read1.".fastq-mcf.fastq";
my $out2 = "";
if($read2 ne ""){
	$out2 = $read2.".fastq-mcf.fastq";
}
my $error_log = $read1.".fastq-mcf.log";

#my $cmd = "$FASTQMCF $ADAPTER -p 10 -k 0 -q 3 -S $read1 -o $out1 ";
my $cmd = "$FASTQMCF $ADAPTER -t 0 -p 10 -l 1 -k 0 -q 0 -x 0 -R $read1 -o $out1 "; 
##the parameter here means, -p 10: Maximum adapter difference percentage
##-t 0: % occurance threshold before adapter clipping. it will always do trimming no matter how small pct of adapter in the file.
##-l 1: All sequence after trimming would be kept, unless it is adapter dimmer
##-k 0: sKew percentage-less-than causing cycle removal, it turned off since it is not good for Bisulfite-seq
##-q 0: do not do quality trimming, since BSMAP would do it, and it will affect our statistics result
##-x 0: 'N' (Bad read) percentage causing cycle removal. it turned off since it will affect our statistics result
##-R: Dont remove Ns from the fronts/ends of reads

if($read2 ne ""){
	$cmd .= "$read2 -o $out2";	
}
$cmd .= "> $error_log \n";
system($cmd);

&extract_result($error_log);
&count_adapter_start_pos($out1);

sub extract_result{
	my $file = shift @_;
	open(FH,"<$file") || die "Fastq-mcf log file $file is missing";
	my $qual_trimmed=0;
	my $adapter_dimmer=0;
	my $adapter_trimmed=0;
	my $total_reads=0;
	while(<FH>){
		chomp;
		my $line=$_;
		if($line=~/^Total reads: (\d+)/){
			$total_reads = $1;
		}
		elsif($line=~/^Too short after clip: (\d+)/){
			$adapter_dimmer += $1;
		}
		elsif($line=~/^Clipped.*Count (\d+)/){
			$adapter_trimmed += $1;
		}
		elsif($line=~/^Trimmed (\d+) reads/){
			$qual_trimmed += $1;
		}
	}
	close(FH);
	if($read2 ne ""){
		$total_reads *= 2;
		$adapter_dimmer *= 2;
	}
	my $qual_trim_perc=sprintf("%.6f",100*$qual_trimmed/$total_reads);
	my $adapter_dimmer_perc=sprintf("%.6f",100*$adapter_dimmer/$total_reads);
	my $adapter_trimmed_perc=sprintf("%.6f",100*$adapter_trimmed/$total_reads);
	print "Total reads: $total_reads\n";
	print "Adapter dimmer: $adapter_dimmer (${adapter_dimmer_perc}%)\n";
	print "Adapter trimmed sequences: $adapter_trimmed (${adapter_trimmed_perc}%)\n";
	#print "Quality score trimmed sequences(base qual <= 2): $qual_trimmed (${qual_trim_perc}%)\n";	
}

sub count_adapter_start_pos{
	my $file=shift @_;
	my $out_startPos = $file.".start_pos.txt";
	my $read_len=50;
	my $linecount_within = 1;
    my $linecount_global = 1;
	open(FH,"<$read1") || die "orignal read file $read1 before trimming is missing";
	while(<FH>){
		chomp;
		my $line=$_;
		if ($linecount_within == 2){
        	 my $read_len=length($line);
        	 last;
        }
        elsif ($linecount_within == 4){
        	$linecount_within = 0;
        }
        $linecount_within++;

	}
	close(FH);
	open(FH,"<$file") || die "read file $file after trimming is missing\n";
	open(OUT,">$out_startPos") || die "can't write on file: $out_startPos\n";
	$linecount_within = 1;
    $linecount_global = 1;
	while(<FH>){
		chomp;
		my $line=$_;
		if ($linecount_within == 1)
        {
            # Double check the format
            print STDERR "Incorrect FASTQ file $file\nLine ${linecount_global}: $line\nMod4 lines should start with \@\n" unless ($line =~ /^\@/);
           
            
        }
        elsif ($linecount_within == 2){
        	 my $start_pos=$read_len-length($line);
        	 print OUT "$start_pos\n";
        }
        elsif ($linecount_within == 4){
        	$linecount_within = 0;
        }
        $linecount_within++;
        $linecount_global++;
	}
	close(OUT);
	close(FH);
}